N.obs <- 5098

# Constructing some simulation data. Instead, you'd use our estimates, 2.5% and 97.5%s bounds, and your own data.
results <- data.frame(lower=runif(N.obs,-3,-1),
                      upper=runif(N.obs,1,3))
results$estimate <- runif(N.obs,min=results$lower,max=results$upper)

# (Check to make sure the simulated estimates are actually within the simulated bounds.)
sum(results$est < results$upper & results$est > results$lower)/N.obs # should = 1, if not run again

# More simulation data.
x <- rnorm(N.obs,0,3)
y <- 2 + 0.25*x - 1.05*results$est + rnorm(N.obs,0,1)
my.data <- data.frame(y=y,x=x,z=results$est,lower=results$lower,upper=results$upper)

# Based on just the estimates
mod1 <- lm(y~x+z,data=my.data)
summary(mod1)

library(truncnorm)
coef.storage <- NULL
N <- 1000 # The larger the N, the better, but it comes at a time cost.
for(i in 1:N){
  # 1. Draw new estimates from a distribution bounded by the lower and upper bounds of the CI.
  new.estimates <- runif(N.obs,min=my.data$lower,max=my.data$upper)
  # 2. Re-run model
  mod.boot <- lm(my.data$y~my.data$x+new.estimates)
  # 3. Store the coefficients
  coef.storage <- rbind(coef.storage,coef(mod.boot))

}

# Estimates are median (or mean) of stored coefficient distribution.
b <- apply(coef.storage,2,median)
lowerCI <- apply(coef.storage,2,quantile,.025)
upperCI <- apply(coef.storage,2,quantile,.975)

# compare
coef(mod1)
b
rbind(lowerCI,upperCI)

# Note that it doesn't get the 'correct' answer, but does correctly identify that it's negative and
# non-zero. (And in most cases, there is no correct answer, since this is simulated data.)